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ABSTRACT 

In a planetary system, a secular particle resonance occurs at a location where the 
precession rate of a test particle (e.g. an asteroid) matches the frequency of one of 
the precessional modes of the planetary system. We investigate the secular interactions 
of a system of mutually inclined planets with a gaseous protostellar disk that may 
contain a secular nodal particle resonance. We determine the normal modes of some 
mutually inclined planet-disk systems. The planets and disk interact gravitationally, 
and the disk is internally subject to the effects of gas pressure, self-gravity, and turbulent 
viscosity. The behavior of the disk at a secular resonance is radically different from that 
of a particle, owing mainly to the effects of gas pressure. The resonance is typically 
broadened by gas pressure to the extent that global effects, including large-scale warps, 
dominate. The standard resonant torque formula is invalid in this regime. Secular 
interactions cause a decay of the inclination at a rate that depends on the disk properties, 
including its mass, turbulent viscosity, and sound speed. For a Jupiter-mass planet 
embedded within a minimum-mass solar nebula having typical parameters, dissipation 
within the disk is sufficient to stabilize the system against tilt growth caused by mean- 
motion resonances. 

Subject headings: accretion, accretion disks — celestial mechanics — hydrodynamics — 
planets and satellites: general — solar system: general — waves 

1. Introduction 

The interaction between a young planetary system and its protoplanetary disk likely plays 
an important role in determining the orbital properties of the planets. Such interactions may be 
important for understanding the observed orbital properties of extra-solar planets (e.g. Marcy 
et al. 1999). Resonances within a gaseous disk likely play a key role in determining planetary 
eccentricities and inclinations. Much of the previous work, starting with Goldreich & Tremaine 
(1980), has emphasized the effects of mean-motion disk resonances, which involve frequencies that 
are comparable to the orbital frequencies of the planets. 
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Another potentially important class of resonances occurs where there is a matching of the pre- 
cession frequency of a test particle (e.g. an asteroid) with the frequency of one of the processional 
modes of the planetary system. The frequencies involved are much lower than in the mean-motion 
case. At such a resonance, the motion of a test particle can be strongly driven by the planets, 
resulting in a high orbital inclination (for a nodal resonance) or eccentricity (for an apsidal reso- 
nance). For example, in the solar system there is an important secular resonance that occurs near 
2 AU due to driving involving Jupiter and Saturn. This resonance is believed to be associated with 
the inner truncation of the asteroid belt (Tisserand 1882). 

It was recognized by Ward, Colombo, & Franklin (1976) that these resonances must have swept 
across portions of the early solar system owing to the effects of the gaseous disk, even if the planets 
do not migrate. The reason is that even a minimum-mass solar nebula can have an important 
influence on the relevant precession rates. As the nebula disperses, the precession rates vary, along 
with the resonance locations (Ward 1981). Regions through which the resonances sweep may be 
driven into significantly eccentric or inclined orbits, such as are observed among the asteroids today 
(e.g. Nagasawa, Tanaka, & Ida 2000). The resonances can also have an important bearing on the 
conditions for planet formation by stirring the planetesimal disk. 

However, collective effects can be important within the disk and may even result in the devel- 
opment of density waves or bending waves. Consequently, the response is less localized than would 
be indicated by test particles, as was recognized by Ward & Hahn (1998) and Tremaine (1998). 
These investigations concentrated on the collective effects of self-gravity in particle disks. 

Secular interactions involving the gas disk of a young planetary system are potentially of 
importance, since much more mass is contained in the gas disk than in the planetesimal disk. 
The disk is a fluid body and is capable of wave motions through (at least) compressive, inertial, 
buoyancy, and self-gravitational forces. It is also capable of dissipating energy through a turbulent 
effective viscosity or through shocks. It is therefore important to understand the dynamical response 
of the disk to secular forcing and the implications for the evolution of the planetary orbits. 

The purpose of this paper is to investigate the outcome of secular interactions in mutually 
inclined planet-disk systems. The relevant waves in the disk are long- wavelength bending waves 
with azimuthal wavenumber m = 1. The equations governing such waves in a protostellar disk are 
fairly well established, while those for long-wavelength eccentric density waves have received less 
attention (although see Ogilvie 2001). The importance of understanding the global m = 1 response 
of the disk to secular perturbations has been emphasized by Tremaine (1998). 

The general outline of the paper is as follows. Sections 2-5 describe the physical model for 
nodal secular interactions and formulate the normal mode analysis. Sections 6 and 7 explore the 
limit of sufficiently slow modes, when the disk responds nearly rigidly. Section 8 describes how 

the effects of mean-motion resonances arc included in the normal mode analysis. Sections 9 and 
10 describe numerical results for Jupiter interacting with the solar nebula. Section 11 provides a 
simple scaling analysis based on a nearly rigid tilt model. Section 12 discusses the relationship of 
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the current approach with that of earher work involving WKB theory. Section 13 describes the 
analysis of a close-orbiting planet in the central hole of a disk. Section 14 discusses Jupiter and 
Saturn interacting with the solar nebula at the viq secular resonance. 

2. Modeling secular interactions 

Planets experience secular gravitational interactions on time-scales much longer than their 
orbital periods. On such time-scales the planets may be considered as continuous rings representing 
their average mass density. These rings are in general elliptical and mutually inclined. Gravitational 
interactions between the rings lead to apsidal and nodal precession of the orbits. 

In the case of small eccentricities and inclinations, the secular evolution of a system of n 
planets can be described in terms of normal modes of a linear dynamical system, n modes for 
eccentricity and n for inclination (e.g. Murray Sz Dermott 1999). The eigenvector of a mode 
describes the relative distribution of eccentricity or inclination among the planetary orbits, while 
the eigenfrequency is the rate at which the pattern precesses. 

In typical situations involving mean-motion resonances (Lindblad or vertical resonances), if a 
particle resonance is located at a certain radius in a fluid disk, the disk responds by launching a 
wave that carries energy and angular momentum away from the resonance. Such waves are usually 
assumed implicitly to be damped through some dissipative process before reaching an edge of the 
disk from which they might otherwise reflect. The gravitational torque between the disk and the 
perturber can then be calculated from a standard formula evaluated at the location of the resonance 
(Goldreich & Tremaine 1979). 

However, in the case of a secular resonance, the frequency may be so small that the wavelength 
would be comparable to (or may even exceed) the size of the disk. In that case the global response 
of the disk must be computed, including explicit dissipation and the correct boundary conditions. 
The exact location of the resonance then ceases to have great importance and the new possibility 
of a global secular resonance arises. 

We focus on the case of inclination, so that the relevant waves in the disk are long-wavelength 
bending waves with azimuthal wavenumber m = 1. As noted above, the equations governing such 
waves in a protostellar disk are fairly well established. 

3. Basic equations 

Let (r, 4>, z) be cylindrical polar coordinates with origin at the central object, of mass M*. We 
consider a protostellar disk of semi-thickness H{r), for which the turbulent viscosity parameter a 
satisfies the condition a < H/r. The linearized equations for bending waves in such a disk have 
been derived in several papers (Papaloizou & Lin 1995; Masset k, Tagger 1996; Demianski k. Ivanov 



-4- 



1997). We present them in the form (Lubow &; Ogilvie 2000) 



dW _ IdG 
dt r dr 



+ T, 



at + H^ T^ + «^^ = ^^• 



(2) 



Here W{r, t) = + ily is the complex tilt variable, which describes the warped shape of the disk. 
Essentially, the disk can be thought of as a continuum of concentric circular rings with radii r and 
unit normal vectors l{r,t). Also g{r,t) = Gx + iGy is the complex internal torque variable. In a 
warped disk there is a horizontal internal torque 27rG(r, t) that is mediated by horizontal motions 
that are proportional to the distance z from the midplane. These stresses are responsible for the 
propagation of bending waves. Similarly T(r, t) = T^+iTy is the complex horizontal external torque 
density acting on the disk. Finally S(r) is the surface density, Cl{r) the orbital angular velocity, 
K{r) the epicyclic frequency and I{r) the second vertical moment of the density. These are defined 

by 

Ti = I pdz, ^— I pz^ dz, (3) 



or 



z=0 



K^=4n^ + 2rn^, 
dr 



(4) 



where $(r, z) is the (axisymmetric component of the) gravitational potential experienced by the 
disk. 



Consider a thin, uniform circular ring of mass nii -C M^, and radius r^, representing either an 
annulus of the disk or the time-average of a planetary orbit. Its contribution to the potential in 
the plane 2; = is 

= -i^ / + rf - 2rn cos <t>r^'^ dct>. (5) 

271" Jo 

Therefore the angular velocity and epicyclic frequency experienced by the ring are obtained from 
summations over all other such rings, in the forms 



2 _ GM^ , 2Gmj 
rf rirj 



Ko{ri,rj) - ^Ki{ri,rj) 



2Gmj 



nr-i 



2r,- 



Ko{ri,rj) ^iCi(rj,r 



(6) 
(7) 



where and K\ are symmetric kernels with dimensions of inverse length, given by 



Km {n , r,- ) = ^ ^ (rf + r] - 2ri cos 4>) "^/^ cos"^ ^ #. 



(8) 



Hereafter we write K\ simply as K. 
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Now consider the tilt interaction between two such rings of masses and rrij, radii and rj, 
and tilt vectors Zj and Ij. For small relative inclinations, the gravitational torque exerted by ring 
j on ring i is 

Tji = GmimjK{ri,rj) k x Ij. (9) 

This expression can be obtained by averaging the torque exerted by a point mass on a circular 
ring (Lubow & Ogilvie 2000) over the orbital motion of the point mass. Similar expressions can be 
found in studies of galactic warps (e.g. Sparke Sz Casertano 1988). In the complex notation, for 
small inclinations from the xy-plane, the total external torque acting on ring i is 

Ti = Y, GmimjKin, rj)i{Wj - Wi). (10) 

The kernels can also be written as 

Km{ri,rj) = ^b^^^ (^^y (11) 

where r> = max(ri,rj), r< = min(ri,rj) and 6^"*^ is the Laplace coefficient. In practice we evaluate 
these in terms of elliptic integrals using Carlson's algorithms (Press et al. 1992). 

These kernels diverge as l^j — r^l ^0. In reality the gravitational interaction remains bounded 
because of the non-zero vertical thickness of the rings, which we have neglected. To take account 
of this in an approximate way, we soften the kernels by replacing 

/ ^ fri+rj-h ri+rj + h\ 
Km[ri,rj) ^ Kml ^ , ^ 1, (12) 

whenever [rj — r^l < h. Here the smoothing length h is to be understood as an approximate measure 
of the thickness of the rings. 



4. Coupled systems of planets and disks 

We consider a general system consisting of multiple planets embedded in a gaseous disk. Pro- 
vided the planets create gaps in the disk, the disk then becomes partitioned into a set of disks. The 
disks and planets interact through gravity. 

We consider planets i = 1, 2, . . . , rip of masses rupi in circular orbits of radii rpj and angular 
velocities ripj. We also consider disks k = 1,2,... ,nd of inner radii and outer radii 6^. The 
inclinations of the planetary orbits are described by the tilt variables Wpi{t). The tilt within a 
disk is described by W^^ (r, t) . We assume throughout that the inclinations are small, so that linear 
theory applies. 

If only secular nodal interactions are considered, the planets can be treated as inclined circular 
rings interacting with each other and with the disks through the torques described in Section 3. 
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We return later (Section 8) to the effect of mean-motion resonances on tilt evolution. However, we 
assume throughout that any evolution of the surface density of the disk, or of the semi-major axes 
of the planetary orbits, may be neglected. 

The secular equations for the planets are then of the form 

"T-pi^pi^pi-^ = ^GmpimpjK{rpi,rpj)i{Wpj -Wpi) 

+ GnipiJ^kKir, rpi)i{Wdk - Wpi) 27rr dr, (13) 

k "•'^ 

while the equations for the disks are 

Efcr^J)^^ = ^^ + Y^GmpiT.kK{r,rpi)i{Wpi-WAk) 

i 

+ r G^k^'^Kir, r')i{W^ - W^k) ^ttt' dr', (14) 
+ ( ] ^Gk + aQQk = — (15) 



dt \ n"^ J 2 " " 4 dr ' 

where = T,i{r') and W^f = Wde{r',t). The final term in equation (14) represents the self- 
gravitation of each disk and the gravitational interactions between disks. 



J2 mpirlPpiWpi + XI / ^kr^^^dk ^irr dr, (16) 

7. J at 



The (complex) total horizontal angular momentum of the system, 

k 

is exactly conserved if, as we assume, the boundary conditions 

Gkiak,t) = gkibk,t) = (17) 

hold at all times t. There is also a conservation law for the vertical angular momentum, which we 
derive in Appendix A. 

5. Normal modes 

We have a set of coupled integro-differential equations that are linear and homogeneous. So- 
lutions may be sought in the form of normal modes, 

Wpi{t) = Wpie'^', Wdk{r,t) = Wdkir)e''^', ^^(r, i) = ^fe(r) e^'^*, (18) 
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where a; is a complex frequency eigenvalue. The real part of w is the precession rate of the mode, 
while the imaginary part is the decay rate. Substituting this and omitting the tildes, we obtain the 
integral and integro-differential equations 

iujmpirliClpiWpi = '^GmpimpjK{rpi,rpj)i{Wpj - Wpi) 

+ J2 f" Gmpil^kK{r, rpi)i{Wdk - Wpi) 27rr dr, (19) 



r dr 



1 

iuj'Ekr'^nWdk = — -r- + ^ GmpiEkK{r, rpi)i{Wpi - W^k) 



+ r G^k^'eK{r, r')i(W^^ - W^k) ^irr' dr' , (20) 

^ J at 



fK^-n^\in^ ikr^n^dWdk 
i^Qk + ^2 ) -^^k + (y-^Qk = (21) 

subject to the boundary conditions 

Qkidk) = Gk{bk) = 0. (22) 

Frequencies O and k are given by equations (6) and (7). All orbits are assumed to be prograde. 
There is always a trivial rigid-tilt mode with a; = 0, = 0, and Wpi = W^^ = constant. 



6. Rigid response 

Secular modes usually have very long periods. If a disk is easily able to maintain radial commu- 
nication (through pressure, viscosity or self-gravitation) on this long time-scale, it will participate 
in such a mode almost as a rigid body, and W will not vary much across the disk. 

Suppose that each disk indeed behaves rigidly (W^k = constant). Then can be eliminated 
by multiplying equation (20) by 27rr, integrating from to b^, and using the boundary conditions. 
We then find 

o^Jp^Wp, = J2 CffiWpj - Wp,) + Cf,\Wak - Wp,), (23) 
^JdkWdk = Yl CfkiWpi - Wdk) + Y CtfiW^e - Wak), (24) 

i ii^k 

where 

is the angular momentum of planet and 

J^k = / Sr^O 27rr dr (26) 
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is the total angular momentum of disk k. The coupling coefficients are defined by 

CfP = GmpimpjK{rpi,rpj), 



^ife^ = / Gmpi^kK{r,rpi)2Trrdr, 

J at 



(27) 
(28) 



/-.dd 



GT,kT,[K{r,r')2TTrdr2iTr' dr', (29) 

Evidently the disks can be treated formally as additional planets if they behave as rigid bodies. By 
introducing Greek indices that run from 1 to rip + n^, wc obtain the algebraic eigenvalue problem 



(30) 



This is formally identical to the multiple-planet problem described by, e.g., Murray & Dermott 
(1999). Since C is a symmetric matrix, we have 



(31) 



Since, further, the coefficients Cp^ are real and positive, the eigenvalues are all real and non-positive. 
Apart from the trivial rigid-tilt mode, all modes precess retrogradely without growth or decay. 

Suppose the system contains only two components, e.g. one planet interacting with a connected 
disk. The rigid solution satisfies 



(32) 



where C = C\2 is the coupling coefficient between the planet and the disk. The eigenvalues and 
eigenvectors are 

a trivial solution corresponding to a rigid tilt, and 



ujJx + C -C ' 




'Wi' 




"0" 


-C UJ2 + C_ 




W2_ 




_0_ 



a; = 0, 



oc 



(33) 



u = —- 



C{ Ji + J2) 
J1J2 ' 



Wi 

W2 



oc 



J2 

-Ji 



(34) 



corresponding to a retrogradely precessing mode in which the tilt of each component is inversely 
proportional to its angular momentum (in order to conserve the total angular momentum). 

In a three-component system, e.g. a planet interacting with interior and exterior disks, the 
rigid solutions satisfy 



UlJl + C12 + Ci3 —C12 — Cl3 

— C12 IJ0J2 + C12 + C23 ^(^23 

— C13 — C23 W J3 -|- C13 -|- C23 





-Wi' 




"0" 




W2 













.0. 



(35) 
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where subscript 1 refers to the planet, and subscripts 2 and 3 to the inner and outer disks. 



In the general case, the eigenvalues and eigenfunctions are algebraically complicated. However, 
consider the case in which the planet is much less massive than either disk. The first non-trivial 
mode is the equivalent of equation (34), but involving the two dominant components: 



UJ ^ —- 



J2J3 



W2 



(X 



-J2 



The planet's tilt is then driven according to 



C12W2 + C13W3 



(36) 



(37) 



W Jl + C12 + Ci3 

Secular resonance occurs here if the frequency is such that the denominator vanishes. 

Secular resonance of this type could also occur, for example, in a system of two massive planets 
and one low-mass disk. In this case the disk would be responding globally as a rigid body. This is 
quite different from the situation in which a secular particle resonance is located somewhere within 
a disk. 



The second non-trivial mode has 

{C12 + C13) 



UJ ^ —- 



Ji 



\w2im\ < \wi\. 



(38) 



Here only the planet's orbit is significantly tilted, and precesses at a rate determined by the two 
disks. 



7. Nearly rigid response 

The internal torque Qt required for a disk to respond rigidly can be determined from equation 
(20). In general, equation (21) will then indicate that a slight warping is in fact required. The rigid 
solution can be understood as the leading term in an expansion of the solution in powers of a small 
parameter. If e = H/r is the angular semi-thickness of the disk, and Cs ~ HVl the isothermal sound 
speed, the small parameter here is a characteristic value of ioj/eVi)'^, or (Lor/cs)"^. This parameter 
is small when the disk can communicate effectively, i.e. when the time taken for a bending wave 
to propagate across the disk is less than the wave period. 

Now suppose the disk behaves nearly rigidly, so that we may pose an expansion in powers of 
this parameter, of the form 

UJ = UJ 

(0) + ^(1) _^ ^ (39) 

Wpi = + + (40) 

W^k = W^f+W^'^{r) + ---, (41) 

Gk = gi^\r) + g^^\r) + .--. (42) 
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At leading order we obtain the rigid solution described in Section 6, denoted by the superscript 
(0). At the next order we find 



+ J2 GmpiEkK{r, r^MW^]^ - W^]^) 27rr dr, (43) 

1 ^r(^) 

i 

+ E r GT.kT.',K{ry)i{W'}i^ -W^^)2T:r' dr', 

n J at 



.(1)^ 

Ak J 



02 



(1) 



subject to the boundary conditions 



ai'^(a,)=ar(M=o. 



(1)/ 



(44) 
(45) 

(46) 



After some manipulations we may eliminate q'^^ by integration and obtain 



a J a 



k "^"fe 



-LO 



(0) 



k2 _Q2\ q 



gf^|2 27rrdr. (47) 



Since does not depend on a but is proportional to S, this shows that nearly rigid modes 
are damped at a rate Im(a;(^)) proportional to a/e^ (for a disk and planetary system of given 
dimensions). The imaginary part of equation (47) is a special limit of the conservation law derived 
in Appendix A. 

It is important to note that the self-gravitational term in equation (20) has no effect on the 
results we have derived. The reason for this is that, if a disk behaves nearly rigidly, the self- 
precession of the disk is negligibly small even if the mass of the disk exceeds that of the planets. 



8. Mean-motion resonances 

So far we have treated only secular interactions between the planets and disks. However, a 
planet also interacts with a disk through mean-motion resonances (Lindblad resonances and vertical 
resonances). The launching of waves at these locations exerts a torque on the planet. If the planet's 
orbit is inclined with respect to the resonant annulus of the disk, this inclination will evolve in time. 

Borderies, Goldreich, Sz Tremaine (1984, hereafter BGT) calculated the rate of change of 
inclination of a satellite due to Lindblad and vertical resonances. They showed that the inclination 
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typically grows exponentially at a rate that depends on the strengths of the resonances present. 
They assumed implicitly that the launched waves are damped before reaching an edge of the disk 
from which they might otherwise reflect. In the case of a gas disk, such waves can damp by means 
of viscous dissipation, radiative damping (Cassen &: Woolum 1996), or wave channeling (Ogilvie & 
Lubow 1999). For narrow gaseous rings between planets, wave damping may well not take place. 
Instead, reflection from the ring edges may occur, which could greatly reduce the effects of the 
mean-motion resonances. For simplicity, we will also assume that the waves are damped very close 
to the resonances and transfer their angular momentum to the disk there. 

To incorporate the effect of mean-motion resonances, we modify our equations as follows. (For 
notational simplicity we here consider only a single planet and a single disk.) The formulae of BGT 
can be interpreted as giving the horizontal component of the resonant torque between the planet 
and the disk. We use this to add terms to the planet's angular momentum equation, in the form 



2o 



j * 



(48) 



where r^j is the radius of the jth resonance, and s^j is the dimensionless strength of the resonance. 
Similarly, we modify the disk's angular momentum equation to 



dt ^ 

3 



(49) 



As noted above, this assumes that the waves transfer their angular momentum to the disk very 

close to the resonances. Provided that the disk is not very cold, the precise location of the wave 
damping is unimportant because the torque is then communicated through the disk over a large 
radial extent. 



The radii and strengths of the resonances are as follows (BGT); here we have neglected the 
shifts of the resonant radii caused by precession. Inner vertical resonance (m > 2): 

2/3 



X = 



m — 1 



m + 1 



Inner Lindblad resonance (m > 2): 



X = — 



m — 1 



m 



2/3 



Outer vertical resonance (m > 2): 



X 



m — 1 



Outer Lindblad resonance (m > 1): 

2/3 



X = — 



Sr = 



24(m - 1) 



rmrx 



12(m - 1) 



2/3 



TTX 



24(m + 1) 



l^ix) 



htjlix) 



n 2 



mTT 



12(m + l) 



(50) 



(51) 



(52) 



(53) 
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These quantities are listed in Tables 1 and 2 for all resonances satisfying |rr — rp| > O.lrp. Note 
that the inner Lindblad resonances act in the opposite sense to the other resonances. 

Although equation (49) formally requires the disk tilt W to have a cusp at each resonance, 
these cusps will be very weak in a disk with good radial communication. For a disk that responds 
nearly rigidly, the effect of the mean-motion resonances is to add a small imaginary part to the 
planet-disk coupling coefficient C^^: 

Gm? 

CPd ^ C-Pd + 27ri -J^Srj^rjTrj. (54) 
j * 

The sum over resonances is always positive for an exterior disk but usually negative for an interior 
disk. According to equation (31), positive sums of this kind tend to cause growth of the non-trivial 
normal modes. This must compete, however, with the viscous decay implied by equation (47). 



9. Model for disk and numerical method 

We adopt a simplified model for the solar nebula. The disk model is similar to that described in 
Lubow & Ogilvie (2000). It is characterized by an angular semi-thickness H/r = e = constant, and 
a surface density S = So/r~^/^, where Sq is a constant and /(r) a function that is close to unity 
in most of the disk, but tapers linearly to zero at the edges (including gaps around planets). The 
width of the tapers is equal to the local value of H. The vertical structure is that of a polytrope of 
index n, which gives I = 'EH'^/{2n + 3) (as defined in eq. [3]) if H is the true semi-thickness of the 
polytrope. The basic parameters of the model are then e, a, and n, together with the dimensions 
and mass of the disk. We normalize the surface density by quoting a nominal disk mass, 

/•rout 

"idisk = / Sor~^/^ 27rr dr = 47^I:or^(^^ , (55) 
Jo 

where rout is the outermost radius. The true disk mass is slightly less than mdisk owing to the 
various edges. 

We consider the nebula in the presence of a Jupiter-like planet. As standard parameters, we 
adopt nip = O.OOIM*, mdisk = O.OIM*, e = 0.075, a = 0.005, and n = 3/2. The planet creates 
a gap in the nebula, leading to two separate disks whose parameters are taken to be ai = O.lrp, 
6i = (1 — g)rp, 02 = (1 + g)rp, 62 = 20rp. The standard relative half-width of the gap is taken to 
he g = 0.2. This situation may be compared with Jupiter in a solar nebula extending out to about 
100 AU. The effect of varying the important parameters will be considered below. 

Numerical solutions to the equations of Section 5 are obtained by discretizing the disk, con- 
verting the integrals into sums, and converting derivatives into centered finite differences, with 
Wdfc defined on a set of rir rings for each disk, and Qk on the boundaries between neighboring 
rings. Planet-disk and disk-disk interactions are computed with softened kernels K^, as defined 
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by equation (12). The smoothing length h = m.ax.{Hi, H2) is taken to be the maximum of the 
semi-thicknesses of the two rings concerned. Equation (49) is implemented by identifying the ring 
in which each resonance falls, and representing the delta function as l/5r, where Sr is the width of 
the ring. The problem then reduces to a generalized eigenvalue problem involving a non-Hermitian 
matrix of dimension np+nd{2n,. — 1). This is solved numerically using the Fortran Sun Performance 
Library routine ZGEGV on a Sun Ultra 10, using rir = 200. Numerical results are quoted in units 
such that G = = rp = 1. The modes are normalized such that Wp = 1. 

10. Numerical results for Jupiter 

We first consider a disk with standard parameters, except that it extends only to 62 = 2 instead 
of 62 = 20. [Accordingly, we set mjisk = (0.1)"'^/^0.01 in this case to fix the same normalization of 
surface density as in the standard model. The true mass of the disk is only 0.00179.] This case is 
easier to understand because the disk is able to maintain good radial communication. 

Only two of the non-trivial modes are of interest. The first mode ('mode 1') has a precession 
frequency of —7.02 x 10"*^. The inner disk has a nearly rigid tilt W « —20 and the outer disk 
W -1-20. If the mean-motion resonances are neglected, this mode has a viscous decay rate of 
4.0 X 10~^; when they are included, the mode acquires a net growth rate of 1.0 x 10~^. [This mode 
corresponds to equation (36) in rigid disk tilt theory.] 

The second mode ('mode 2') has a precession frequency of —1.46 x 10~^. The inner disk has 
a nearly rigid tilt W ^ —0.64 and the outer disk W « —0.75. If the mean-motion resonances are 
neglected, this mode has a viscous decay rate of 7.1 x 10~^; when they are included, the mode 
acquires a net growth rate of 6.6 x 10~^. 

The remaining non-trivial modes all have significantly larger precession rates and damping 
rates (with or without mean-motion resonances). They are the proper bending modes of the disks, 
modified by gravitational coupling to the other components of the system. 

When we proceed to the standard model, the numerical solution becomes significantly more 
complicated. The outer disk extending to 62 = 20 is not able to maintain good radial communication 
for typical precession rates. The precession frequencies estimated from the nearly rigid theory are 
comparable to the frequencies of global bending modes in the outer disk. As a result, the nearly 
rigid theory is no longer accurate for the outer disk. 

To simplify the presentation of results, we focus on three modes that are of greatest interest 
in that they involve the least warping, and have the smallest precession frequencies and damping 
rates. 'Mode 1' and 'mode 2' arc related to the nearly rigid modes mentioned above, while 'mode 
3' derives from a bending mode of the outer disk that interferes with the other two when the outer 
disk fails to behave rigidly. The structures of the three modes are plotted in Fig. 1. The modes 
include the effects of mean-motion resonances. 
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In Figs 2 and 3 we plot the precession rates and decay rates of these three modes for the 

standard model, as we vary each of the parameters a, e, m^isk^ 9 about its standard value. 
Notice that the variation of decay rate with some parameters is non-monotonic. In particular, the 
decay rate of mode 1 peaks at mdisk ~ 0.02. 

In Appendix A, we derive an expression for the local decay rate of a warped disk due to 
dissipation. In Fig. 4 we plot for mode 1 the decay rate, as given by equation (A3). Notice that the 
outer disk dissipation dominates over the inner disk dissipation. Although the dissipation peaks 
near the planet, it is broadly distributed throughout the outer disk. In the outer disk, the warp 
\rdW/dr\ is greatest near the disk radial midpoint and is broadly distributed. Quantity \rdW/dr\ 
must vanish at the inner and outer disk edges, due to boundary condition (22). The reason that the 
peak in Fig. 4 in the outer disk occurs near the planet is that the magnitudes of various base-state 
quantities, such as S and 0,, are largest there. 

In summary, an inclined planet gives rise to large-scale warps. The damping effects of secular 
interactions dominate over tilt excitation by mean-motion resonances in our model system, provided 
a > 0.001. For the standard model, the decay time of the longest lived mode, mode 1, is about 
2.7 X 10^ yr. Although this time-scale is short compared to typical disk lifetimes of several million 
years, it is comparable to or longer than the expected planetary migration time-scale at this radius 
(Lin, Bodenheimer, & Richardson 1996; Ward 1997). 



11. Simple estimates 

We compare the numerical results with the nearly rigid theory of Section 7 for the standard 
disk model. To do this, we carry out some rough estimates and drop all factors of order unity. 
Consistent with our standard model, we assume below that the mass of the planet is less than, or 
comparable to, the masses of the two disks, and also that the angular momentum of the outer disk 
is greater than, or comparable to, the angular momenta of the planet and inner disk. 

The coupling coefficients are dominated by the parts of the disk closest to the planet, where 
|r — Tpl/rp ~ g. Here the kernel may be estimated as K ^ \/{g^r). Thus 

CP<^ ~ ^GmpSr (56) 

and 

~ GSV^. (57) 

Here S is a typical surface density near the planet. The precession rate of a mode can be estimated 
from either equation (36) or equation (38) as 



(58) 
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We then estimate G from equation (20) as 

g ~ iio^r^nW. (59) 
We apply these estimates to equation (47) to obtain an estimate of the viscous damping rate, 

When mean-motion resonances are taken into account, the planet-disk coupling coefficient 
acquires a small imaginary part, as in equation (54). The sum over resonances can be estimated 
using arguments adapted from BGT. For g <^ 1, the number of resonances in the disk scales as 
m ~ 1/g. The strengths of the Lindblad resonances, which dominate here, scale as s-c ^ m^. Thus 

. 1 Gmillr 

^From equation (36) this would provide a resonant growth rate 

M-)~^=5f^a (62) 

The ratio of viscous decay to resonant growth can then be estimated as 

decay f f 9^ 



growth 



The first factor depends on the viscosity and temperature of the disk. The second is related to the 
ratio of the disk mass within (say) two gap widths of the planet to the mass of the planet itself. 
For our standard parameters, these first factor is close to unity, while the second is somewhat less 
than unity. This favors tilt growth, provided that the disk indeed behaves nearly rigidly. 



12. Connection with WKB bending wave theory 

Disks subject to forcing by a misaligned companion have been studied in the context of plan- 
etary rings, but with strikingly different results. Within Saturn's A ring, tightly wrapped waves 
are launched from a mean-motion vertical resonance by a misaligned satellite, as described by Shu, 
Cuzzi, Sz Lissauer (1983). Furthermore, the torque carried by the waves is independent of a and e, 
unlike the present case. 

As emphasized earlier, the difference lies in the fact that the waves here are of low frequency 
such that the parameter ur/cs is typically less than unity for protostellar disks with secular reso- 
nances. For the mean-motion vertical resonances in planetary rings, this quantity is much greater 
than unity. 
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Nonetheless, a formal connection between the equations used in this paper and those used for 
planetary rings can be made by considering an artificially cold protostellar disk, such that cjr/cg 
is much greater than unity. Furthermore, we consider the case that the disk mass is much greater 
than the planetary mass. Appendix B outlines the derivation of the dynamical equations in the 
WKB limit. In that limit we obtain 

where cr^ = ±1 with the sign depending on whether the waves are trailing or leading, and 

2Gm 

D{r) = 2wO + —^K{rp, r). (65) 
For a; O, we can re-express D{r) in the standard form 

D{r) = -{n' -uf, (66) 

where 

,2 92$' 



or 

z=0 



(67) 

z=0 



Potential is due to the star and axisymmetric contributions of the planet, but does not include 
contributions due to the disk. This equation is identical in form to equation (20) of Shu et al. (1983) 
for m = 1, apart from the pressure term, the third term in the equation (64). 

The pressure term arises because the warping of the disk introduces horizontal pressure gradi- 
ents. The resulting shearing horizontal motions give rise to the hydrodynamic torque 2ttG, which 
is the dominant factor in the propagation of bending waves in protostellar disks. The warping is 
described by the n = mode in the notation Lubow and Pringle (1993). The dispersion relation 
for adiabatic oscillations of this mode in an exactly Keplerian disk, whose unperturbed state is 
vertically isothermal, follows from equation (54) of Lubow and Pringle (1993) with dimensionless 
frequencies F = (fi — uj)/Q ~ 1 and k = 1. In the low-frequency regime (|a;| <C 0,), one obtains 
then that u = zizkHQ/2, independent of the exponent for adiabatic compression 7. Although the 
form of the dispersion relation is reminiscent of that for a compressive acoustic mode, this mode is 
largely incompressible, as suggested by the lack of dependence on 7 in the dispersion relation. The 
same dispersion relation (as was obtained by Papaloizou & Lin 1995) follows from equation (64), 
with W varying as exp {ikr), and I = {cs/Cl)'^T,, and by taking G = (ignoring gravity), k = J7, 
and a = 0. 

For this limit where the planet mass is small compared to the disk mass, the resonance condition 
D{r) = becomes 

where {K) denotes the mass weighted average of K over the disk (see Appendix B for more details). 
For the standard disk model, the resonance condition is satisfied at radial distances of about 0.7 
and 1.3 within the inner and outer disks respectively. 
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This system contains a secular particle resonance where the precession rate of a particle matches 
the precession rate of the disk-planct system. A free particle (e.g. an asteroid) will precess due to 
the gravitational effects of both the disk and the planet. That is, a free particle resonance will occur 
when D{r) = in equation (66), but ^' now has an additional contribution due to the disk. The 
particle resonance location is therefore not the same as the disk resonance location. In other words, 
in locating a particle secular resonance, the gravity of the gas disk must be included. Particles 
oscillating vertically in a non-oscillating gas disk have an additional vertical restoring force that 
alters their precession frequency and shifts their resonance site. This is a special property of the 
m = I bending modes (tilt modes), which is a consequence of the fact that there is no self-precession 
in a rigidly tilted disk. Near resonance, a tilt mode essentially behaves rigidly, due to its locally 
long wavelength. 

In the case of a planetesimal (particle) disk, the nature of the pressure and viscous interaction 
may well differ from that of a gaseous disk. On the other hand, the secular resonance condition 
given by equation (68) likely applies to a low-mass planet interacting with a cold planetesimal disk. 

13. Planet in central hole 

Current models for the early evolution of close-orbiting extra-solar planets, such as 51 Peg b 
(Mayor & Queloz 1995) with an orbital radius of about 0.05 AU, involve the planet migrating into 
a central hole in the disk (Lin et al. 1996). The planet is envisioned to migrate into the hole by 
disk tidal forces until its 2:1 outer Lindblad resonance lies just inside the inner edge of the disk. 
The planet would then remain at this radius, without further migration. 

We consider the tilt stability involving both secular interactions and mean-motion resonances 
of such a planet-disk system, based on the standard model parameters in Section 9. In this case, 
there is only an outer disk of mass 0.01 M*, which extends from just outside the planet's 2:1 outer 
Lindblad resonance l.Slrp (0.066 AU) to outer radius 2000 rp (100 AU). Other disk parameters are 
the same as in the standard model (e = 0.075, a = 0.005, and n = 3/2). The only mean-motion 
resonance present in the disk is the rather weak 1:3 resonance (2.08 rp), which acts to increase 
inclination. 

Fig. 5 shows the eigenfunction of the lowest mode of the system. The eigenfrequency is —5.11 X 
10"^ + 4.6 X 10~^i in units of Qp, corresponding to a decay time-scale of only 3.9 x IC^ yr. The 
inner part of the disk interacts nearly rigidly with the planet, owing to good radial communication 
oJVp/cs -C 1 (see lower panel of Fig. 5). However, the outer part of the disk is not in good 
communication [cor /cg is of order unity) , resulting in a warp and dissipation. 
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14. Disk and two planets 

We consider here the case of Jupiter and Saturn interacting with the portion of the solar 
nebula that lies interior to the orbit of Jupiter, using the equations of Section 5. This case is of 
interest because the classical uiq nodal secular resonance lies within the disk. We ignore the effects 
of mean-motion resonances. 

The rigid tilt model involving three objects (here the two planets plus disk) was considered at 
the end of Section 6. As is mentioned there, a secular resonance of the two planets with the disk is 

possible, if the disk is low in mass compared with the planets. If the disk mass is greater than the 
planet masses, then the planets individually interact with the disk in a manner similar to the case 
in Section 10. 

We have computed the modes of the system numerically. We consider Jupiter (planet 1) and 
Saturn (planet 2) in circular orbits with the present values of mass (mpi = 0.0009545, mp2 = 
0.0002858) and semi-major axis (rp2 = 1.833). A test particle would experience nodal secular 
resonances at semi-major axes 0.381 and 2.39, corresponding to 1.98 AU and 12.4 AU (Murray 
&; Dermott 1999). We consider a partially depleted solar nebula interior to Jupiter's orbit. Our 
standard model has mdisk = 0.001, e = 0.075, a = 0.005, n = 3/2, and disk inner and outer radii 
ai = 0.1 and bi = 1 — g = 0.8, respectively, in units of Jupiter's orbital radius. The inner secular 
resonance for a test particle therefore lies inside the disk. 

For these parameters, the numerically determined eigenfrequencies of the two lowest modes 
are —2.27 x 10""^ -|- 6.0 x 10~^i and —7.66 x 10~"^ + 5.1 x 10~''i. The response of the disk is very 
nearly rigid and nothing special happens at the location of the secular particle resonance. This is 
to be expected, since wr/cg ^ 1 throughout the disk. The effect of self-gravitation is negligible 
even though the mass of the disk is comparable to that of the planets. 

If a low-mass disk is extremely cold, it is unable to establish radial communication very ef- 
fectively on secular time-scales. It cannot respond nearly rigidly, but instead launches a wave at 
the location of the particle resonance. (The rigid response is essentially the limit of a wavelike 
response when the wavelength becomes large compared to the size of the disk.) In Fig. 6 we show 
the wavelike response generated at the i^iq secular resonance at 2 AU. A very thin disk (e = 10~^) 
is required to see this effect, and even then the wavelength is not much shorter than the size of the 
disk (so a WKB or tight- winding treatment of the wave might not be particularly accurate). We 
set a = e = 10~^ so that the wave would be partially damped on reaching the inner boundary. We 
also chose a low disk mass, mjisk = 10^^, so that the disk has a Toomre parameter Q > 1. The 
eigenfrequency for this lowest mode is —2.35 x 10^"^ + 1.6 x lO^^i, corresponding to a damping time 
of about 1.1 X 10^ yr. The resonantly launched waves damp inclination faster than those in the 
warmer disk model, even though the disk mass is much lower. 
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15. Summary 

We have carried out a normal mode analysis of planet-disk systems that includes the effects of 
secular interactions and mean-motion resonances. The planets and disks interact gravitationally, 
and the disks communicate internally through gas pressure, self-gravity, and turbulent viscosity. 
Secular interactions of misaligned planct-disk systems give rise to global effects in the disks, owing 
mainly to the effects of gas pressure. The low frequencies u) associated with secular modes allow 
pressure to communicate effectively over radial distances in disks that are comparable to (and 
somewhat greater than) the planet's orbital radius. Over such distances, disks behave rigidly with 
little dissipation. On the other hand, over larger distances (~ Cs/oo), which may arise in a continuous 
disk exterior to the planet's orbit, the radial communication breaks down and large-scale warps 
occur, along with enhanced dissipation (see Figs 1 and 5). The protostcllar disk response to secular 
interactions typically falls between that of a global nearly rigidly tilted disk (see Section 7) and 
that of a locally launched wave (see Sections 12 and 14). 

Secular interactions act to decrease inclination by means of turbulent dissipation of horizontal 
shearing motions within the warp for a simple a disk prescription. The dissipation of the warp is 
peaked somewhat near the planet, but is broadly distributed throughout the disk (Fig 4). Secular 
interactions are in competition with the effects of mean-motion resonances (BGT), which can act 
to increase inclination. For standard disk parameters, the secular interactions involving a single 
planet generally dominate and the inclination decays in time. For a young Jupiter interacting with 
the solar nebula, tilt decay takes places for a > 0.001, with a typical alignment time-scale of order 
10^ yr (see Fig. 3). For a close-orbiting planet in the central hole of the disk, the alignment is more 
rapid. These results suggest that a planet formed within a disk will remain coplanar with the disk, 
at least as a consequence of disk interactions. 

The classical viq particle secular resonance, due to Jupiter and Saturn, is expected to have 
resided within the solar nebula. Again, due to pressure effects, this resonance is greatly broadened 
in the gas disk. Consequently, the nebula responds nearly rigidly, with some decay in inclination. 
A very cold planetesimal disk might exhibit a genuinely wavelike response (cf. Fig 6). 

The analysis in this paper has focused mainly on the effects of a single planet interacting with 
its interior and exterior disks. For a system of a multiple planets that opens multiple disk gaps, 
the situation is more complicated because of reflections that occur at the disk edges associated 
with the gaps. Such gaps can decrease the effective size of the disk region over which a planet can 
interact, thereby resulting in a weaker warp and weaker dissipation. The hydrodynamic torque 
cannot propagate across a gap, while the gravitational torque may at least partially do so. It is 
quite possible that much longer lived noncoplanar modes may result if the outer part of the nebula 
is almost disconnected in this way. 

This work was supported by NASA grants NAG5-4310 and NAG5-10732, the STScI visitor 
program, and the Institute of Astronomy visitor program. GIO acknowledges the support of Clare 
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A. Derivation of the local decay rate of bending disturbances 

Starting from equations (13)-(15) we obtain the conservation law 



where 



k •'"fc 

The quantity is the vertical angular momentum associated with the bending disturbance. When 
an initially horizontal ring of matter having orbital angular momentum J is tilted through an angle 
(3 = \W\, the change in the vertical component of angular momentum is J(cos /3 — 1) ~ — ^ JjW^I^ 
when \W\ ^ 1. The term proportional to represents an additional angular momentum stored 
in the bending wave. For a short-wavelength (WKB) bending wave in a non-self-gravitating disk, 
both terms in the integral for —Lz are equal. For a nearly rigidly tilted disc, the first term in 
the integral dominates. In the appropriate limit, this first term agrees with the expression for the 
angular momentum of a bending wave given by Bertin & Mark (1980, eq. [CI 2]). 

In an inviscid system is conserved. When viscosity is present, the positive definite quantity 
—Lz decays monotonically to zero. The negative angular momentum of the bending disturbance is 

transferred to the disk through viscosity, and causes accretion. (It is nevertheless consistent that we 
neglect the time-dependence of S in a linear theory, because the induced accretion rate is quadratic 
in W.) 

We now define a local damping rate 7(r) by 

^ k 

This has the property that 

In the case of a normal mode, the radial average of 7 corresponds to the damping rate Im(a;) of 
the mode. 



7 
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B. Derivation of the bending wave equations in the WKB limit 

In this Appendix we derive equations (64) and (68) of the text. We consider equations (19)-(21) 
for a single planet interacting with a single disk. 

iujmprlnpWp = J Gmp^K{r,rp)i{W - Wp)2Trr dr, (Bl) 

iuo^r^nW = -^ + GmpJ:K(r,rp)i(Wp-W) 
r dr 

+ J GT.T.'K{r, r')i{W' - W) 27rr' dr', (B2) 



We now apply the WKB approximation for a cold disk, wr S> Cg. In this regime we assume 

d\nW 



dr 



> r-^. (B4) 



The requirement on ^ is used in the derivation of the warp equations (see Papaloizou & Lin 
1995; Lubow and Ogilvie 2000). Prom equation (B3), we obtain 

r ~ ~4 [w + (k2 - 02)7(20) - iaa\ dr^ ' ^ ' 

For the self-gravity term, we have that 



/ 



GSS'ii:(r, r')i{W' - W) 27rr' dr' « iTrGSV^^^. (B6) 



To obtain equation (B6), wc recognize that self-gravity is dominated by local contributions in the 
WKB limit and apply the approximation that K{r,r') ~ r/[27r(r' — r)^] for r' ~ r. The integral is 
computed by extending r' to the complex plane and integrating over a closed contour. The contour 
begins along the negative real axis, includes a small semicircle around the singularity at r' = r, 
continues along the real axis, and closes on itself through a large semicircle. The large semicircle 
resides in either the upper or lower half-plane, depending on whether the waves arc trailing or 
leading. Equation (B6) then follows from the residue theorem. This choice of trailing or leading 
waves determines the sign of the result. Substituting equations (B5) and (B6) into equation (B2), 
we obtain equation (64). 

Integrating equation (B2) over the disk area and using equation (Bl), we have that 

Ja{W) = -JpWp, (B7) 



where {W) denotes the angular momentum weighted average of W. In the limit that Jp <C Jd, we 
expect that W on the right hand side of equation (Bl) can be ignored after integration. It then 
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follows from equation (Bl) that 



jG'EK{r,rp)2Trrdr _ 



u = = - — . (B8) 



'p^'p "^p 



Applying this equation for co to the resonance condition D{r) = 0, for D defined by equation (65), 
we obtain the resonance condition (68). 
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Table 1. Mean-motion resonances interior to the planet 
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Table 2. Mean-motion resonances exterior to the planet 
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Fig. 1. — Eigenfunctions of the tilt variable W for the three lowest modes, with mode 1 plotted at 
the top, mode 2 in the middle, and mode 3 at the bottom. The plots are normalized such that the 
the tilt of the planet (dot) is unity at a radius of unity. The solid curves are for the real part of the 
eigenf unction, while the dashed curves denote the imaginary part. For the particular disk model 
adopted (with constant -ff/r), W{r) is also equal to the vertical displacement from the midplane in 
units of the local disk scale- height, if the planet resides at one disk scale-height above the midplane. 
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Fig. 2. — Variation of the precession rates of the three lowest modes as the parameters a, e, mdisk) 
and g are varied independently about their standard values. The solid, dashed, and dotted lines 
correspond to modes 1, 2, and 3 referred to in the text. 
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Fig. 3. — Variation of the decay rates. The plot follows the notation of Fig. 2. 




Fig. 5. — Eigenfunction of the tilt variable W for the lowest order mode involving a close-orbiting 
planet in the central hole of a disk. The plots follow the notation of Fig. 1. 
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Fig. 6. — Eigenfunction of a secular mode for an artificially cold disk with a = 10~^, e = 10~^, and 
^disk = 10~^. Real (solid) and imaginary (dashed) parts of the complex tilt variable arc shown. 
The mode is normalized such that VFjupiter = 1, while Wsaturn ~ —2.47. A wave is launched at the 
location of the particle secular resonance (r = 0.381) due to the Jupiter-Saturn interaction. 



